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ABSTRACT 

A theoretical engineering thermal analysis of man-made "permafrost 
islands" as offshore drill sites is conducted by developing a two- 
dimensional numerical model from the heat conduction equation. The 
model is formulated by solving an equivalent variational principle 
based on the finite element method. Phase transformation due to 
freezing and thawing is assumed to occur isothermally to promote 
solution economy. Solutions in the time domain are obtained by 
the Crank-Nicolson method. 

Boundary conditions may be arbitrarily imposed or simulated at the 
air-surface interface. Surface boundary conditions are simulated by 
using a one-dimensional surface energy exchange model which is coupled 
to the two-dimensional model. Local meteorological data is employed 
for the surface energy simulation. Successful tests were conducted 
that compared the numerical model solutions with exact solutions to 
freezing front propagation. Also, comparisons were made between model 
results and in situ measurements of soil temperatures at Fairbanks, 
Alaska. Finally, the model was applied to the problem of offshore 
"permafrost island" drill sites. 
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Chapter 1 
INTRODUCTION 

1»1 Alternatives to Conventional Offshore Drilling Structures 

Recent studies indicate that offshore regions in the Beaufort and 
Chukchi Seas may contain extensive oil and gas reserves. Some geologists 
believe that offshore areas tend to be more productive than adjacent 
land areas. If this is true for arctic waters, the implications are 
staggering considering the recent petroleum discoveries along the arctic 
coast. It can be expected that exploratory offshore drilling in the 
Beaufort Sea will commence sometime within the forseeable future with 
initial probes restricted close to shore. 

Offshore drilling experience in hostile arctic environments has 
been limited to Cook Inlet, Alaska. Although annual growth of sea ice 
in Cook Inlet is only 3 feet, design wind, wave, and earthquake forces 
on offshore structures have been found to be negligible when compared 
with expected lateral ice forces. Obviously, the coastal environment 
in the Beaufort Sea is much more hazardous to ocean structures since 
annual sea ice reaches thicknesses of 7 to 9 feet and is subject to 
continuous motion due to the driving forces of wind and current. 

Other dangers include the possibility of pressure ridges, ice islands, 
or multi year ice striking a structure. 

Assuming an adequate knowledge of critical design loading due to 
ice, the application of known engineering techniques for offshore drill- 
ing structures could probably be physically implemented but, at an 
economically prohibitive expense. The arctic thermal regime may permit 
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consideration of two additional alternatives for development of offshore 
drill sites: 

1. Construction of man-made permafrost islands from dredged sea- 
floor materials that are frozen into an icy matrix strong 
enough to resist lateral forces exerted by annual sea ice 
movement. 

2. Construction of structurally reinforced man-made ice islands 
grounded at the desired location. 

1.2 Scope and Objective 

This report is addressed to the first alternative. The primary 
problem, thermal analysis of nearshore man-made permafrost islands, is 
considered. The objective is to develop a versatile engineering tool 
based on known thermodynamic theory and make appropriate assumptions to 
facilitate engineering application. A mathematical model is formulated 
for computer solution utilizing the finite element method. Numerical 
solutions are compared with analytical solutions prior to application 
to example permafrost island problems. This report is restricted to a 
theoretical study relying on available data for soil thermal properties 
and available historical records for meteorological conditions. The 
results of several long-term simulations are presented. They are 
considered to be of value to engineers desiring design information on 
equilibrium freezing front propagation in various soil structures 
subjected to offshore arctic meteorological conditions. 
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Chapter 2 

FINITE ELEMENT FORMULATION OF TRANSIENT HEAT CONDUCTION WITH PHASE CHANGE 

2.1 Freezing and Thawing of Soils 

Heat transfer with phase change due to freezing or thawing is the 
principal thermal process of interest in applied arctic soil engineering 
problems. When attempting to mathematically quantify or simulate the 
physics of these processes, the engineer must assess what level of model 
sophistication is necessary in the context of the proposed application. 
The primary factors affecting his decision are the extent to which the 
physical properties of the soil system are known and the quality of 
information available for determining external boundary conditions. 

As a general rule, the physical parameters of the soil system will 
dictate how the phase change problem should be solved. Three alterna- 
tive methods are considered for modeling the latent heat of fusion 
problem for saturated soils: 

A. Phenomenono logical 

B. Thermodynamic 

C. Analytical - isothermal 

Each method and its inherent limitations is briefly discussed. 

2.1.1 Phenomenonological Modeling 

It is well known (Keune and Hoekstra, 1967) that certain moist or 
saturated soils do not undergo isothermal phase change. Lukianov and 
Golovko (1957) proposed a phenomenonological equation to describe such 
a process in a soil-water-ice system: 
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3 3cf> 3<j> 

3^0® = Pgt CD 

where x = distance from the air-soil interface (ft) 
t = time (day) 

4> = temperature C°F) 

K = K (x,t,<J>) = thermal conductivity (BTU/ft-day-°F) 

F = apparent volumetric heat capacity (BTU/ft^-°F) = C - L dG/d<j> 

Cjq = (x,t,<j>) = volumetric heat capacity of the soil -water-ice 

system 

G(<J>) = Ice content (lb/ft^) 

L = latent heat of fusion (BTU/lb) 

Equation (1) mathematically describes the freezing-thawing phenome- 
non in a one-dimensional zone of finite thickness. The freezing of a 
small differential element is accompanied by a temperature dependent 
release of latent heat which decreases the time rate of change of 
internal energy if net heat extraction from the element remains invariant. 
The phenomenonological approach is equivalent to introducing an "apparent 
volumetric heat capacity" as defined above. 

Although phenomenonological modeling of phase change remains 
analytically intractable, it is possible to numerically solve equation (1). 
Nakano and Brown (1972) developed a one -dimensional finite difference 
scheme from equation (1) and satisfactorily simulated heat conduction 
in tundra soils near Barrow, Alaska. It is important to stress that 
the success of any phenomenonological model is heavily dependent on 
the following factors: 
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1. Unfrozen soil water content as a function of temperature must 
be known or experimentally determined. 

2. Numerical simulation of phase change requires small time incre- 
ments to avoid solution instabilities caused by the injection 
or removal of latent heats. 

3. Applied boundary conditions must be known on a micro-scale. 

2.1.2 Thermodynamic Modeling 

Thermodynamic modeling of soil phase change also incorporates 
equation (1) to predict soil column temperatures as a function of time 
but does not rely on experimentally determined soil freezing curves to 
simulate the temperature dependent evolution of latent heat. Recent 
research (Koopmans and Miller, 1966) has show that many soil moisture 
characteristic curves and soil freezing curves are similar. Since a 
soil moisture characteristic curve defines the relationship between 
water content and soil pore water pressure, detailed calorimetric 
analyses to determine unfrozen water content in soils as a function of 
temperature can be avoided if temperature-pore water pressure relation- 
ships during phase change are known. 

It is possible to calculate soil water pore pressure in soils under- 
going phase transition by assuming that pure ice and water can only 
co-exist in thermodynamic equilibrium as defined by the Clausius - 
Clapeyron equation for a bulk fluid (Sears, 1964) 

LdT 

dp= T( v "-v') ( 2 ) 

where L = latent heat of fusion 

T = existing absolute temperature 



v" = specific volume of ice 

v' = specific volume of water 

dp = incremental change in pressure 
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It is sometimes assumed that the Clausius - Clapeyron equation applies 
to unsaturated soils and that ice pressures in the soil pores are at 
atmospheric pressure. It is probable that this equation does not apply 
to saturated soils without modification. 

A typical numerical procedure for thermodynamic simulation of soil 
freezing is to compute a soil temperature profile from equation (1), use 
this to obtain soil water pore pressure by numerically integrating 
equation (2), and enter the soil moisture characteristics curve to 
determine the change in water content. The net change in water content 
over a given time increment determines the magnitude of latent heat 
that must be released in the following time step. A new or updated 
apparent volumetric heat capacity is computed to continue the temporal 
solution of equation (1). This method of discretely modeling phase 
transformation requires the utilization of small time increments to 
avoid unstable oscillations in computed temperatures. Furthermore, 
soil moisture characteristic curves must be experimentally determined. 

2.1.3 Analytical Modeling of Isothermal Phase Change 

Analytical solutions commonly employed in the engineering analysis 
of phase change problems are simplified versions of the theory of ice 
growth introduced by Franz Neumann in the 1860's. Examples are Stefan's 
(1889) and Berggren's (1943) applications of the Neumann theory to ice 
growth and frost penetration problems. Since any theoretical formulation 
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of phase change (commonly referred to as the "problem of Stefan") is 
nonlinear, classical solution techniques are not applicable (Carslaw 
and Jaegar, 1959). 

Because of the mathematical complexities involved, the practicing 
engineer is usually forced to seek a published analytical (exact) 
solution that represents the best approximation to his physical problem. 

A significant disadvantage of mathematically exact solutions is the 
limitation imposed by simplifying assumptions and solution space geo- 
metry. To analytically solve the isothermal phase change problem, it 
is usually necessary to assume a semi-infinite isotropic medium initially 
at a uniform temperature (i.e,. , simplified boundary conditions are 
chosen so that the governing equation of state can be solved) . The 
imposed boundary condition is a sudden change in temperature that subse- 
quently induces phase change in the medium. It is the restrictive 
imposition of simple boundary conditions that severely limits the 
general application of exact solutions. 

Several Russian investigators have been successful in analytically 
modeling more realistic boundary conditions. Kolesnikov (1946) presented 
a derivation which incorporated the major hydrologic energy exchanges 
affecting the growth of sea ice, and more recently, Portnov (1962) 
developed an exact solution for isothermal phase change with time 
dependent boundary temperatures. 

Because of the relative paucity of exact solutions and the complexi- 
ty of modern engineering problems, most research efforts have been 
directed toward achieving approximate methods of solution. The most 
effective have been finite difference and finite element methods which 
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rely on analytical solutions as a means of testing numerical techniques 
prior to simulating physical problems that remain analytically insoluble. 
Hence the primary importance of known theoretical solutions is that of 
investigating the accuracy of a numerical procedure whose validity has 
not been theoretically established. 

2.2 Derivation of a Two-Dimensional Finite Element Model 

An approximate numerical solution to the parabolic differential 
equation of transient heat conduction in two dimensions is developed by 
employing the finite element technique based on the so-called Rayleigh - 
Ritz formulation. Phase change is assumed to occur isothermally to 
avoid previously discussed difficulties common to thermodynamic and 
phenomenonological methods of solving the latent heat problem. 

2.2.1 Variational Finite Element Formulation 

The variational functional approach for the finite element formula- 
tion of transient heat conduction relies principally upon the theoretical 
contributions of Biot (1970) who introduced a unified variational 
analysis of heat transfer analogous to the energy theorems of classical 
mechanics. Given the partial differential equation of two-dimensional 
transient heat conduction in Cartesian coordinates 




( 3 ) 



where <j> = temperature 

K x and fC, are t 
directions. 



.nd 1C, are the thermal conductivities in the x and y 
directions, respectively 



G = G(x,y,t), a variable heat source or sink 
pc = volumetric heat capacity 
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t = time 

the transition to an equivalent variational statement is sought by- 
applying an extension of the fundamental lemma of variational calculus 
(Weinstock, 1952). The theorem is applied by noting that if in a given 
time interval, the functional 

x = //f(x,y,*,|l|i)dxdy 

R 9x 9y (4) 

is to be minimized over a bounded region R, a necessary and sufficient 
condition is that at any given instant of time the unknown function 
4>(x,y,t) satisfy the Euler equation 



3 3f 
"5x^9 (9<J>/9x) 



0 



+ I r . 

9y l 9(94>/9y) J 




(5) 



subject to <j> obeying the same boundary conditions of equation (3). The 
problem is to obtain a functional which when substituted into equation 
(5) reduces to equation (3) . 

Assuming remains invariant, an equivalent variational formula- 
9t 

tion of equation (3) is obtained if the functional 

x s iV {KxC ir )2 + + 2(pc lt ■ G ^> dxd y (6) 

can be minimized for all admissible states of <J>. For equation (6) to 
be a true minimizing principle, the second variation of the functional 
with respect to the unknown state variable must be greater than zero. It 
can be shown that this is indeed true for transient heat conduction 
(Myers, 1971). It is less restrictive to require the functional to be 
stationary (i. e. , only the first variation is zero) and demonstrate 



that the numerical solution will converge to an exact solution. 
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Equation (6) is formally solved by discretizing the domain R into 
numerous subregions or "elements" and assuming that the temperature 
distribution within each element can be approximated by an n^ order 
polynomial. In this report, a linear polynomial is selected to 
approximate the state variable in a two-dimensional finite element mesh 
composed of arbitrarily shaped triangular elements. 

Upon substitution of a linear temperature distribution shape 
function into equation (6), the first variation of the functional is 
evaluated by applying Leibnitz's rule for differentiation of integrals. 
Thermal properties are assumed to be constant within each element although 
abrupt changes in properties between elements are permitted. Summation 
of the integrated contributions of all elements produces a system of 
simultaneous linear first order differential equations of the form 

[K]{<t>) + [C]{|£} = (7) 

where [K] = system conduction matrix 
[C] = system capacitance matrix 
{G} = system load column matrix 
{<J>} = nodal temperature column matrix 
The system conduction and capacitance matrices are banded and symmetric, 
significantly reducing computer memory requirements for storing arrays. 

2.2.2 Incorporation of Isothermal Phase Change 

Simulation of heat conduction and phase transformation in two 
dimensions involves a minor modification of the variational formulation 
of transient heat conduction. As previously discussed, the solution 
domain is discretized by a finite element mesh composed of triangular 
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elements. The mesh geometry is stored in computer memory by reading in 
nodal points, their respective global coordinates, thermal properties 
of each element (frozen and unfrozen), and moisture content of each 
element. After this has been accomplished, numerical simulation of 
phase change may commence. - 

Consider a common node point "A" in a segment of a two-dimensional 
finite element mesh (Figure 1) whose temperature is dependent upon any 
arbitrary boundary condition imposed along S. Phase change is not 
permitted to occur at node "A" until a requisite amount of latent heat 
has been exhausted or added. The prevailing quantity of latent heat 
for node "A" is determined by adding one third the weight of material 
(subject to phase transformation) of each contributing element to a 
nodal latent heat accumulator array assigned to node "A". 

As an example, consider the cooling and freezing of node "A". 

After each time increment that a solution is being generated, the temp- 
erature of node "A" is scanned to determine whether or not it has 
dropped below a prescribed point where phase transformation is assumed 
to occur (TFAZ) . If the node has not been previously frozen and the 
computed temperature has dropped below TFAZ, then the quantity of 
latent heat evolved during the previously computed time step becomes: 
DELQ = C U (TFAZ - Computed Temperature) * Volume (8) 

where = a weighted volumetric heat capacity based on the heat 

capacities of the soil-water mixture assigned to node "A M 
that must undergo phase change. 

Volume = volume of soil-water assigned to node "A" 

The quantity DELQ is subtracted from the latent heat accumulator for 
node "A" and if more latent heat must still be exhausted, the computed 
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Fig. I 



Typical Finite Element Mesh 
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temperature at the node is brought back to TFAZ prior to obtaining a 
new solution at the end of the next time increment. This is analogous 
to simulating phase change as an isothermal process. The time rate of 
latent heat generation is governed solely by the solution of the heat 
conduction equation and remains insensitive to large time steps. 

If the required amount or more than the required amount of latent 
heat has been extracted, node "A" and the volume of soil-water assigned 
to it are considered to be frozen. The thermal properties of the 
elements common to node "A" are then updated based on the volumetric 
proportions of soil -water and soil -ice present in the respective 
elements. If an excess amount of latent heat is removed, the residual 
is used to calculate how far below TFAZ the temperature at node "A" 
should be. A weighted Cf is used for this computation. All weighted 
volumetric heat capacities (Cu and Cf) are computed from the thermal 
properties of elements common to nodal points undergoing phase change. 

The thawing process is handled in a similar manner. Provisions are 
made to monitor which nodes are frozen or thawed so as to avoid re- 
freezing frozen nodes. Boundary nodes remain unaffected by the checking 
routines that scan nodal temperatures. After phase transformation at 
a node occurs, the latent heat accumulator assigned to that node is 
recomputed. Recomputation permits the simulation of long-term cyclic 
freezing and thawing. 

2.2.3 Imposition of Boundary Conditions 

The variational functional expressed by equation (6) automatically 
incorporates boundary conditions of either specified temperature or 
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zero heat flux. Inclusion of additional modes of heat transfer that 
may occur over various portions of the boundary surface, shown in 
Figure (1), is accomplished by adding appropriate surface integrals 
to equation (6) which upon minimization yield the desired boundary 
conditions. The surface integrals for specified heat flux, convection, 
and radiation heat transfer across the boundary are 

S q<Ms * if K(-t> 2 - 2<M)ds * 0/ (ie B * 5 " (9) 

B f 2 B c Br 5 

where q = specified heat flux across the boundary 
<{> = boundary node temperature 
<p m - ambient temperature near the boundary 
K = convective heat transfer coefficient 
a = Stefan-Boltzmann constant 
eg = Emissivity of the boundary surface 
e^ = Emissivity of the ambient surroundings 
Outgoing heat fluxes are considered positive. 

The combination of equations (6) and (9) and their subsequent 
minimization produce a generalized equivalent variational statement 
describing transient heat conduction with its associated boundary 
conditions. Substantiating mathematical details are provided in 
Appendix A. Minimization of surface boundary integrals produces terms 
that are combined with the system conduction and load matrices of 
equation (7) . Nonlinear radiation heat transfer terms give rise to 
solution difficulties since unknown temperatures are introduced into the 
system matrices. Methods of circumventing this problem include lineari- 
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zation or modeling net radiation exchange as a specified heat flux. 

These methods are discussed in greater detail in section 3.2.2. 

2.2.4 Solution in the Time Domain 

Linearization or removal of nonlinear terms from the system 
matrices allows preservation of the previously discussed system of 
ordinary differential equations represented by 

[K]{<f>} + [C]{||} = {G} (10) 

and greatly simplifies the time domain solution. After a considerable 
amount of experimentation with various numerical time domain solution 
techniques, the Crank-Nicolson method was selected since it provided 
the most satisfactory results in terms of optimizing time step sizes, 
of providing acceptable solution accuracy, and of minimizing numerically 
induced oscillations. 

The Crank-Nicolson method utilizes the arithmetic mean of the 
derivative of the state variable at the beginning and end of each time 
step to move the solution ahead in time. Thus 

(<j>} v+1 = 4) v + ♦ (4>} v+1 ) (11) 

where At = time step size 

v = a reference point in time 

Upon premultiplication by the capacitance matrix and substitution of 
equation (10), equation (11) becomes 

CM + ^[C])4> V+1 = C^[C] - [K]){<J»} V + 2{G} (12) 

Since (<f>} V is the initial or previously computed temperature 
distribution and the system matrices [K], [C] , and {G} remain constant 
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in each time step, equation (12) can be further simplified to the 
matrix equation 

[W]{<J>} V+1 = {Z> (13) 

where [W] = a symmetric matrix in band form. 

Instead of employing a time consuming matrix inversion algorithm, the 
solutions, {<}>} v+1 , are computed by the Gaussian elimination method which 
exploits the symmetric and banded structure of the [IV] matrix. 

2.3 Comparison of Finite Element Solutions with Analytical Solutions 

Involving Phase Change 

Because the method formulated to solve the latent heat problem is 
a numerical approximation, comparison of finite element solutions with 
several exact solutions to simple isothermal phase change problems is 
justifiable. Two problems are considered: one in which thermal para- 

meters remain unaffected by phase change and secondly, one where phase 
transformation produces a significant variation between frozen and 
unfrozen thermal properties. 

2.3.1 Simulation of Isothermal Phase Change with No Change in 
Thermal Parameters 

As a first example, the problem of the solidification of a homo- 
geneous slab of finite thickness is considered. An exact solution 
discussed by Allen and Severn (1952) is used for comparison purposes. 

In Figure 2, a two-dimensional finite element mesh for analyzing the 
problem is shown. Figures 3 and 4 summarize the finite element results 
obtained from two different time domain solutions: implicit and Crank - 
Nicolson. Good agreement between both numerical solutions and the exact 
solution is evident. Note that the Crank-Nicolson time domain solution. 
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Figure 2. Solidification of a slab of liquid: finite element 

mesh. 

Distance between nodes 1. and 19 = L = 60. AX = L/6 = 10 

K = 1; pc = /2 (thermal properties remain constant) 

Uniform initial temperature = 45° 

Temperature of phase transformation = 45° 

Latent heat of transformation = 70.26 pc/unit volume 

The entire perimeter with the exception X = L is 
thermally insulated. 

The end X = L at t = 0 + is dropped to and subse- 
quently maintained at a temperature of 0° . 



Figures 3 and 4. Solidification of a slab of liquid: progress of 

freezing front through the slab. 

Inset diagram illustrates nomenclature in each 
figure. 

The freezing front is denoted by line F at a typical 
time t, and the frozen zone is shown shaded. 

Exact solution for the equation of the solidification 
locus is: 

(X-L) 2 = p|(1.07)t 

as obtained from Jones' theory (Allen and Severn, 
1952) 



/ 
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Fig. 2 

Solidification of o Slob of Liquid; Finite Element Mesh 
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Fig. 3 

Soiidif icotion of o Slab of Liquid : Progress of Freezing Front through Slob. 
Implicit Method Utilized for the Time Domoin Solution. 
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Solidificotion of o Slob of Liquid: Progress of Freezing Front through Slob. 

CronK-Nicolson Method Used for the Time Domoin Solution 
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Figure 4, permits similar solution accuracy with time steps twice as 
'large as those employed in the implicit time domain solution. Figure 3. 

This example demonstrates the accuracy of finite element discretiza- 
tion of the freezing process. Judicious selection of time steps and 
finer mesh subdivision will improve solution accuracy with increased 
computer simulation time. A much more rigorous and realistic test is 
to consider the following theoretical problem where substantial changes 
in thermal properties of the conducting medium occur upon phase change. 

2.3.2 Simulation of the Neumann Theory of Freezing Front Propaga- 
tion in a Saturated Soil Column 

The application of the Neumann theory to the prediction of frost 
penetration depth in soils is well known. A transcendental equation is 
solved to obtain a constant of proportionality whose magnitude is 
dependent upon the thermal properties (frozen and unfrozen) of a semi- 
infinite soil column, initial uniform temperature distribution, and 
imposed initial constant boundary temperature. This constant yields the 
exact solution for freezing front propagation when multiplied by the 
square root of elapsed time. 

Comparison of a one -dimensional finite element simulation with a 

Neumann solution for freezing front propagation is shown in Figure 5. 

The Neumann solution was extracted from an example cited by Jumikis 

(1966). The freezing front locus is plotted versus the square root of 

time to demonstrate the stability of the finite element solution and 

enable computation of the numerical constant of proportionality. The 

theoretical constant of proportionality is .575 ft/day - ^/^ and the finite 

- 1/2 

element solution constant of proportionality is .561 ft-day" ' . With 





V 
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Figure 5. 



Comparison of a variational and analytical solution 
to the Neumann problem. 

A. Physical constants : saturated sand at a void 

ratio of .50; specific gravity = 2.65. 

B. Thermal parameters : 

1 . Unfrozen : 

Kq = 25.7 BTU/ft-°F-day 
C u = 42.7 BTU/ft£°F 

2. Frozen : 

Kf = 32.2 BTU/ft-°F-day 
Cf = 29.3 BTU/ft -°F 

C. Initial and imposed boundary conditions: 



Initial uniform temperature distribution of 
36°F. A constant boundary temperature of 14°F 
is maintained. Isothermal phase transformation 
is assumed to occur at 32°F. 



D. Solutions : 

Analytical and variational solutions are: 
x = .575 /t (theoretical) 
x = .561 /t (finite element) 
where x = frost penetration depth (feet) 
t = time (days) 

Variational solution was obtained by employing 
a one -dimensional finite element mesh. Spacing 
between nodes was one foot with numerical 
computations conducted in .5 day increments. 



23 



L 

e 




o 



Comparison of o Variatianai Solution and on Analytical Solution to the Neumann Problem. 
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such a small discrepancy in constants of proportionality, it would 
take an elapsed time of 15 years before the error between computed and 
theoretical frost penetration depth exceeded one foot. At that time the 
theoretical freezing front would have penetrated 42.6 feet below the 
ground surface. 
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Chapter 3 

SIMULATING AIR-GROUND INTERFACE ENERGY EXCHANGES 

3.1 Methods of Approximating Heat Transfer at the Air-Ground Interface 

The thermal processes occurring at the interface between a soil 
system and the atmosphere are often too complex for an exact descrip- 
tion; yet, their very existence is crucial to conducting an analysis 
of the thermal regime. The dominating processes are absorbed incoming 
solar radiation, net long wave radiation exchange between the atmosphere 
and surface boundary, turbulent heat transfer, evaporation, snow and 
vegetative cover effects, and the soil temperature gradient at the air- 
ground interface. There are three alternative methods that may be 
employed to simulate the net effect of these processes: 

1. Impose a soil surface boundary temperature based on known 
ambient air temperature or in^ situ temperature measurements 
of the soil surface. 

2. Compute an equivalent temperature acting across a known 
thermal resistance to produce a net flux into the soil surface. 
With this method soil surface temperatures are computed rather 
than imposed. 

3. Mathematically model each of the individual thermal processes 
as precisely as possible and solve the resulting system of 
nonlinear equations simultaneously with equation (7). 

The method formulated in this report for simulating soil surface 
energy exchanges is a hybrid combination of alternatives 2 and 3. 
Individual thermal processes are modeled on a macroscopic scale, con- 
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verted to an equivalent temperature, and coupled to a variable turbulent 
heat transfer coefficient to simulate heat fluxes transmitted across the 
air-surface interface. This approach generates surface temperatures 
and eliminates the necessity of solving a system of nonlinear equations. 

3.2 Modeling the Individual Processes of Energy Exchange 

From a macroscopic viewpoint an energy balance at the air-ground 
interface can be approximated by the expression 

Qsoil + Qswr + Qnlwr + Qconv " Qevap = ® 
where Q S oil e heat flux into the soil column 

Q swr = short wave radiation component 
Qnlwr = net l° n S wa ve radiation exchange term 
Qconv = iavhulent heat transfer component 
Qevap = evaporative heat loss term 
and where all heat fluxes are expressed in BTU/ft -day. 

Prior to discussing how equation (14) is incorporated into the variation- 
al analog for transient heat conduction, the methods employed to deter- 
mine each of the above energy components are briefly discussed. 

3.2.1 Short Wave Radiation 

Incident short wave radiation is determined by 

1. Applying CRREL nomographs (Gerdel et.ad. , 1954) for computing 
available solar energy as a function of time and latitude. 

2. Modifying available solar energy with empirical relationships 
(Eagleson, 1970) which take into account altitude of cloud 
base and per cent cloud cover. 
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3. Selecting appropriate surface albedos from Wechsler and Glaser 
(1966). 

3.2.2 Long Wave Radiation Exchange 

As discussed in section 2.2.3, the net radiant heat emission term 
is nonlinear. If the surface emissivities of two radiating bodies are 
the same, a linearization technique discussed by Mettler (1966) and 
Pivovarov (1973) could be applied to provide 

Q n lwr = 4 a£4>£(4>-<U (15) 

where all terms are as defined in section 2.2.3. 

Since the emissivity of a soil or snow surface is different than 
that of the atmosphere, it is easier to compute net long wave radiation 
exchange from the Stefan-Boltzmann law and apply net radiation heat 
transfer as a discrete flux across the boundary surface. An "apparent 
atmospheric emissivity" is assigned from empirical formulas (Eagleson, 
1970) which incorporate per cent cloud cover, cloud base altitude, and 
atmospheric water vapor pressure as influencing variables. 

Inaccuracies are introduced because previously computed boundary 
surface temperatures are used to determine net radiative heat transfer 
in the next time increment of the transient heat conduction problem. 

This is preferable to formulating an exact variational statement because 
computer solution by iteration is avoided. 

3.2.3 Evaporative Heat Losses 

Evaporative heat transfer is not directly simulated but is estimated 
from data and applied as an outgoing flux across the surface boundary. 

The magnitude of evaporative flux for a given region is determined by 
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analyzing the potential evapotranspiration (PET) or actual evapotranspira- 
tion (AET) losses from data given by Patric and Black (1968). A time 
dependent evaporative heat flux function is constructed with a temporal 
shape closely approximating that of available solar energy. Its inte- 
grated area equals the selected PET or AET. In this report, it is 
desired to simulate the effect of evaporation from a saturated soil 
system; hence, PET is employed to compute an evaporative heat flux 
function. Evaporation during the winter months is assumed to be neglig- 
ible. 



3.2.4 Turbulent Heat Transfer 

Turbulent heat transfer is modeled according to the relationship 
^conv “ ^ ^surface " ^air^ 

where K = a turbulent heat transfer coefficient that is a 

function of wind velocity, surface roughness, and 
atmospheric stability 

^surface = surface temperature 

<}> a i r * ambient air temperature 

The theoretical and experimental contributions of Scott (1964, with 
subsequent verification in 1969) are used to compute a turbulent heat 
transfer coefficient that is allowed to vary at the end of each time 
increment or when meteorological parameters change. 

3.2.5 Snow Cover Effects 

It is not desired to simulate the hydrologic morphology of a^snow 
cover over the soil surface but to observe the gross insulating effect 
on a seasonal basis. In this respect heat conduction through an average 
winter snow cover is modeled using the previously discussed energy 
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exchanges. Snow cover density is estimated from data presented by 
Biello (1969) and thermal conductivity from curves discussed by Mellor 
(1964). 

Minor sophistications are included to permit degradation of snow 
surface albedo when the ambient air temperature rises above 32°F. 
Simulation of heat conduction through the snow cover is terminated 
when the computed snow surface temperature is greater than or equal to 
32°F. 



3.3 Finite Element Simulation of Surface Energy Exchange 

Having determined individual heat flux components, net energy 

exchange across a two-dimensional surface is simulated by coupling 

one-dimensional elements having a thermal conductivity equal to the 

turbulent heat transfer coefficient and having a volumetric heat 

capacity of zero to the two-dimensional elements representing the 

boundary surface. An equivalent temperature is calculated from 

<j> = <f> + Sswr. J_9nl ,wr ~ .Qe.vap (17) 

eq amb K 

where = equivalent temperature 

^amb = am ^^ ent a ^- r temperature 

K = turbulent heat transfer coefficient 

and other terms remain as previously defined. The computed temperature 

from equation (17) is then imposed at the end of each one-dimensional 

element. 



/ 

The imposition of equivalent temperatures at discrete nodes produces 
a net uniform flux (satisfying equation 14) across the two-dimensional 
boundary surface once they have been properly inserted into the system 
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matrix equations. Although substantiating mathematical details are 
omitted for the sake of brevity, note that if the equivalent tempera- 
ture is equal to the ambient air temperature, only turbulent heat 
transfer is simulated. Figure (6) i s a schematic presentation of the 
concept . 

3.4 Simulation of Surface Energy Exchange under Fairbanks Field Conditions 

Numerical simulations of surface energy exchange were conducted 
for Fairbanks, Alaska summer and winter climatic conditions to permit 
evaluation of several solution techniques. Numerical instabilities 
encountered in modeling winter snow cover effects led to the development 
of the equivalent temperature approach for handling thermal processes 
at the air-surface interface. 

3.4.1 Surface Energy Simulation: 1-10 July 1972 

Using local meteorological data for Fairbanks, summer surface 
temperatures of exposed silt were simulated over a ten day period. Silt 
thermal properties were computed from the experimental data of Kerstens 
(1949). A clear day outgoing evaporative heat flux was determined from 
Fairbanks evapotranspiration loss data (Patric and Black, 1968) and 
further modified by a factor equal to the ratio of actual solar energy 
received to theoretical available clear day solar energy. Input data 
and results are shown in Table 1. Computed ground surface temperatures 
were approximately 8°F warmer than the ambient air temperature. 

t 

3.4.2 Surface Energy Simulation: 20-29 January 1973 

A ten day winter simulation of heat conduction in a soil column 
with an 18 inch snow cover was conducted to allow comparison of 
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SYMBOLS 

<£e<j = Equivalent Temperature 

= One - Dimensional Surface 
Energy Exchange Element 

k = Turbulent Heat Transfer Coefficient 

pc = Volumetric Heat Capacity=0 



Fig. 6 

Modeling Surface Energy Exchanges 
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I. Meteorologicol Doto 



Day 


Clear Day 
Solar Flu* * 


Clear Day 
Evap. Flu* 


Avt. Daily 
Temperature(°F) 


Relative 
Humidity (%) 


Ave. wind Speed 
(HPH) 


Cloud Base 
(1000 Feel ) 


Claud Cover 
(Tenths) 


1 


2140. 


750. 


62. 


71. 


6.8 


34. 


8. 


2 


2140. 


750. 


62. 


59. 


5.3 


34. 


4. 


3 


2140. 


750. 


! 67. 


00 

• 

• 


4.0 


34. 


10. 


4 


2140. 


750. 


70. 


• 

00 


7.1 


34. 


9. 


5 


2140. 


750. 


68. 


55. 


6.3 


21.6 


7. 


6 


2121. 


770. 


69. 


53. 


5.3 


34. 


4. 


7 


2121. 


770. 


72. 


38. 


9.9 


34. 


3. 


8 


2121. 


770. 


70. 


35. 


8.1 


34. 


1. 


• 9 


2121. 


' 770. 


70. 


41. 


7.5 


34. 


3. 


10 


2121. 


770. 


69. 


51. 


5.3 


34. 


1. 



I. Simulation Results 



Doy 


Simulated 
Transmitted 
Solar Radiation 


Simulated 
Evaporative 
Heat Loss 


Simulated Net 
Long Wave 
Energy Lass 


Simulated 
Turbulent 
Heat Loss 


Simulated 
Sail Surface 
Temperature (°F) 


1 


1813. 


635. 


139. 


500. 


70. 


2 


1816. 


636. 


331. 


400. 


71. 


3 


1811. 


634. 


169. 


415. 


79. 


4 


1812. 


634. 


136. 


553. 


79. 


5 


1439. 


502. 


134. 


424. 


79. 


6 


moo: 


653. 


322. 


398. 


78. 


7 


1800. 


653. 


301. 


476. 


78. 


8 


1802. 


654. 


458. 


390. 


75. 


9 


1801. 


654. 


• 

VO 

8 


437. 


77. j 


10 


1802. 


654. 


453. 


367. 


77. 



X All Energy Terms ore Expressed in BTU /FT 2 - Day 



Toble I 

Summer Simulation ( I - 10 July 1972) 
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measured and computed temperature profiles at the end of the 10 day 
period. Of secondary interest were the computed snow surface tempera- 
tures. Average Fairbanks January snow density and thermal properties 
were determined from values presented by Wheeler Q973). In this 
simulation, solar and long wave radiation fluxes were placed directly 
across the snow surface without employing an equivalent temperature 
and turbulent heat transfer coefficient. Solution instabilities caused 
by day to day fluctuations in meteorological parameters were encountered 
when all surface energy balance terms were simulated in such a manner. 
Available options to solve the stability problem were those of utiliz- 
ing very small simulation time increments or applying the equivalent 
temperature concept. The latter was chosen. 

Results of the winter simulation are shown in Table 2. Ten 
elements were used to simulate the snow cover and six elements to model 
a 45 inch soil column. A lower soil boundary temperature of 32°F was 
imposed since field observations indicated a mean January position of 
the groundwater table at a depth of 45 inches below the ground surface. 
Calculated and measured results were in good agreement. 



I. Meteorological Data 
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Day 


Cleor Ooy 
Solor Flu* X 


Clear Ooy 
Evap. Flux 


Ave. Oaily 
Temperature (°F ) 


Relative 
Humidity (%) 


Ave. Wind Speed 
(MPH) 


Cloud Bos* 
(1000 Feel) 


Cloud Cover 
(Tenths) 


1 


93. 


0. 


-36. 


69.7 * 


0.0 


16.3 


9. 


2 


93. 


0. 


-39. 


66.0 


.4 


19.0 


8. 


3 


118. 


0. 


-37. 


68.7 


0.0 


11.3 


9. 


4 


118. 


0. 


-36. 


68.5 


1.3 


11.0 


8. 


5 


137. 


0. 


-38. 


59.5 


.6 


18.7 


8. 


6 


137. 


0. 


-42. 


65.2 


.4 


23.7 


7. 


7 


155. 


0. 


-37. 


68.2 


2.7 


34.0 


3. 


8 


155. 


0. 


-16. 


61.6 


.9 


18.6 • 


7. 


9 


214. 


*0. 


11. 


71.8 


3.5 


15.J2 





10 


214. 


0. 


16. 


73.1 


11.4 


9.0 


9. 


E. S 


C 

o 

a 

D 
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Results 




' » 

HE. Computed and Measured 
Column Temperatures at the 


Day 


Simulated Snow 


Simulated Net Simulated Ground 
Lcr.g Wave (Surface Temp.(°F) 

[Raoiation Exchonpe 




of the 10 Day V/ inter- Simulal 


Surfoce Temp.f 0 ?' 








lul » e 1 1 #e A 


1 


-30. 


77.2 


22.3 




Soit Column 
Depth (Inches) 


S i mu toted 
Soil Temp. (°F! 


Soil Temp. (°F1 


2 


-34. 


132.7 


21.8 




0. 


20.6 


20.5 


3 


-30. 


58.2 


21.5 




2. 


20.8 


21.2 


4 


-30. 


77.1 


21.7 




4. 


21.1 


21.4 


5 


-33. 


127.7 


21.6 


- 


8. 


21.7 


22.6 


6 


-39. 


153.8 


21.0 




16. 


23.3 


24.3 


7 


-44. 


254.3 


19.9 




24. 


25.2 


25.7 


8 


-20. 


129.9 


18.7 




36. 


28.5 


29.5 


9 


5. 


104.3 


18.6 




45. 


32.0 


32.0 


10 


11. 


73.3 


20.6 











* All Energy Terms are Expressed in BTU /FT -Day 



Toble 2 

Winter Simulation (20- 29 January 1973) 
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Chapter 4 

TWO-DIMENSIONAL PERMAFROST ISLAND THERMAL ANALYSIS 
4.1 Statement of the Problem 

A theoretical engineering analysis was conducted of the thermal 
feasibility of man-made permafrost islands. Two different problems 
were considered: 

1. The thermal degradation of an island drill site built on 
location during the winter months and completely frozen 
prior to spring breakup. 

2. The propagation of a winter freeze front through an island 
constructed from dredged seafloor sediments during the summer. 

Simulated island geometry was a right cylinder with a vertical 
height of 25 feet and radius of 24 feet immersed in 20 feet of water. 
Computer storage limitations prohibited investigation of larger island 
geometries. Two-dimensional finite element discretization of a vertical 
cross section of the island is shown in Figure (7). Symmetry in the 
transient heat conduction solution was assumed, thus, only the left 
half of the island is shown. 

Daily meteorological data from Barter Island, Alaska, was utilized 
to simulate surface boundary temperatures. Barter Island is located 
2 miles from the north shore of the Alaskan mainland in the Beaufort 
Sea. The climate is determined by the surrounding ocean environment 
which prevents extremely low winter temperatures and limits the warming 
effects of continuous daylight during the summer months. 
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Two - Oimensionol Finite Element Mesh 
96 Elements 
61 Nodes 
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4.2 First Thermal Analysis: Prefrozen Drill Site 

The first thermal analysis conducted assumed a prefrozen drill 
site built during the winter at an initial uniform temperature equal to 
the ambient air temperature (15°F). Thermal simulation was for the 
period 1 May 1970 through 31 August 1971. 

The upper 4 feet of the island consisted of saturated gravel with 
a porosity of .45 and the remaining 21 feet were saturated silt with a 
porosity of .5. Soil thermal properties were determined from Kersten's 
(1949) report. Gravel in the upper part of the island that thawed 
during the summer was assumed to drain to a moisture content of 10 per 
cent of gravel unit dry weight. Draining would deter thawing, and 
gravel thermal properties were modified accordingly. Subsequent winter 
freezing was assumed to occur isothermal ly at 100 per cent soil satura- 
tion at 32°F. 

4.2.1 Applied Boundary Conditions 

To conserve use of computer time, upper surface boundary conditions 
were estimated by conducting a one -dimensional finite element surface 
energy exchange simulation. A vertical cross section of the island was 
employed. Barter Island meteorological data was provided in 24 hour 
increments, and daily surface temperatures were computed for use in the 
two-dimensional model. Outgoing summer evaporative heat flux was 
estimated from Figure (8). The presence of a 15 inch mean winter snow 
cover at a density of .33 g/cm was simulated for the period 1 October 
1970 to 15 May 1971. Monthly averages of simulated energy terms are 
shown in Figure (9), where average net flux across the boundary surface 



EVAPORATIVE HEAT FLUX (BTU/FT* DAY) 
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Fig. 8 

Estimated Barter Island, Alaska Evaporative Heat Flux 
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2 Net Flux Across the Surfoce Boundory 

3 Net Long Wave Radiotion Energy Exchonge 

4 Turbulent Heat Transfer 

Fig. 9 

• Simuloted Surface Energy Terms 
15“ Winter Snow Caver Assumed 
Sorter Island, Alosko 
( I Moy 1970 -31 Aug. 1971 ) 
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is derived as a residual term. In Figure (10), simulated average 
1 monthly surface temperatures and snow-ground interface temperatures 
are shown along with average measured ambient air temperatures' for 
comparison purposes. 

The effect of a .02°F per foot geothermal gradient acting across 
a sea bottom of frozen silty clay was simulated and found to be neglig- 
ible. For this reason, the bottom of the two-dimensional island was 
considered to be nonconducting. 

A mean annual seawater temperature of 28.5°F was applied along 
the completely immersed vertical sides of the island. 

4.2.2 Simulation Results 

Two-dimensional computer simulation results are summarized in 
Figures (11) through (14). These figures show computed isotherms for 
the 1970-71 winter and following summer for the cases of mean seasonal 
and negligible winter snow cover. The second case assumes that snow is 
removed to promote maximum soil cooling. In both simulations, approxi- 
mately 2.5 feet of the surface gravel thawed by the end of the summers 
of 1970 and 1971. To determine porosity effects on thaw, a one- 
dimensional simulation predicted a 5 foot depth of thaw if gravel 
porosity was .30 instead of .50. 

4.3 Second Thermal Analysis: Dredged, Unfrozen Drill Site 

A second two-dimensional thermal analysis was undertaken to 
determine winter freezing front propagation in a dredged artificial 
island with the same geometry and soil structure as the previously 
discussed case. Natural extraction of heat was enhanced by assuming 
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Averoge Monthly Simulated Surface Temperotures 
ond Meosured Ambient Air Temperatures 
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Computed Isotherms (°F) on I Moy 1971 
Simulotion Time Span: I Moy 1970 - 31 Aug 1971 

15" Mean Annual Winter Snow Cover Included. 
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Computed Isotherms (°F) on 31 Moy 1971 
Simulation Time Spon : I May 1970 - 31 Aug 1971 
15" Mean Annual Winter Snow Cover Included. 
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Computed Isotherms (°F) on I May 1971 
Simulation Time Span: I May 1970 - 31 Aug 1971 
Winter Snow Cover Assumed to be Negligible. 
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Computed Isotherms (®F) on 31 Aug. 1971 
Simulation Time Span: I Moy 1970 - 31 Aug 1971 

Winter Snow Cover Assumed to be Negligible. 



46 



that no significant amount of snow would be allowed to accumulate on 
'the surface of the island. 

Salinity of the seawater in the saturated silt and gravel -was 
assumed to be 15 parts per thousand which, depressed the freezing point 
to 30.4°F. Oceanographic measurements (Hufford, 1974) indicate that 
salinities may very along the arctic coast from 3 to 20 parts per 
thousand in the upper 50 feet of water. 

4.3.1 Imposed Boundary Conditions 

As before, a one-dimensional thermal analysis of surface energy 
exchange was accomplished utilizing a saturated silt soil column to 
determine if enough silt would freeze over one winter to justify a two- 
dimensional thermal simulation. Results are shown in Figure (15) where 
depth of freezing front penetration is plotted against degree days of 
frost. It was judged that enough silt did freeze to consider a two- 
dimensional analysis worthwhile. 

It was also recognized that the imposition of a mean annual sea 
temperature of 28.5°F would not be completely representative of imposed 
winter boundary conditions along the vertical sides of the island. For 
this reason, the growth of sea ice was simulated and the temperature 
distribution in the ice was imposed along the upper submerged portion 
of the island considered to be in contact with the ice. Geothermal heat 
flux effects were again assumed to be negligible. 

4.3.2 Simulation Results 

Two-dimensional simulation results are shown in Figure (16). 
Approximately the upper 13 feet of the island froze during the first 
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Fig. IS 

One - Dimensional Simulation of the Freezing 
of Saturated Silt Using 
Barter Island Cimatologicol Data 
(15 SEP. 1970 - 15 May 1971) 
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Two - Dimensional Simulation of the Freezing 
of o Soturoted Grovel and Silt Island 
Computed Isotherms (°F) on 15 Moy 1971 
Simulation Time Span : I Nov. 1970 - 15 May 1971 
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winter. To assist in predicting the elapsed time for the entire 
^island to freeze, a one-dimensional simulation incorporating geothermal 
heat flux was run with the results shown in Figure 0-7) • It was 
estimated that the island would be completely frozen after 5.5 years. 
Extrapolation of results indicated an equilibrium frost penetration 
depth of 27 feet would be reached in a saturated silt soil structure 
with a 4 foot saturated gravel over-burden. 

A second one-dimensional simulation using saturated gravel pre- 
dicted that a gravel island of the same dimensions would be completely 
frozen in 3.5 years. Computed equilibrium freeze front penetration in 
gravel subjected to Barter Island climatological conditions was 48 feet 
after an elapsed time of 18 years. Simulation results are shown in 
Figure (18). 

4.4 Discussion of Results 

Both thermal simulations demonstrated that surface degradation of 
island drill sites due to summer thawing would be less than 3 feet, 
which is quite minimal. The primary explanation for this is that 
average Barter Island summer cloud cover is approximately eight-tenths, 
and consequently the magnitude of transmitted short wave radiation is 
low. As an example, average simulated transmitted solar heat flux, 
from Figure (9) , during the month of June was only about 50 per cent 
of the available theoretical clear day solar energy. 

The convergence of computed isotherms CFig ures (.11) and (13)) near 
the air-water interface indicated that a large quantity of heat was 
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Fig. 17 

Simulated Elapsed Time to Reach an Equilibrium 
Freeze Front in a Saturated Gravel /Silt 

Soil Structure 
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Fig. 18 

Simuloted Freeze Front Propogotion- in Soturoted Grovel 
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extracted from the surrounding water, transmitted through the island, 
and exhausted to the atmosphere. during the winter. 

Of considerable concern is the intrusion of heat along the vertical 
face of the island during the summer months as shown by the position of 
the 29°F isotherm at the end of the 1971 summer ([Figure (14)). The 
danger of thawing can be eliminated by constructing prefrozen island 
drill sites with low salinity water. Dredged islands would not present 
a major problem since brine cells should migrate towards the freezing 
front during the winter and be rejected into the unfrozen silt. The 
high salinity silt would then become diluted during the summer months. 
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Chapter 5 

CONCLUSIONS AND RECOMMENDATIONS 



5.1 Conclusions 

The objective of this report was to formulate and develop a viable 
numerical model designed to solve arctic engineering problems concerning 
transient heat conduction accompanied by isothermal phase change. The 
first phase of the work involved comparison of the numerical model with 
several analytical solutions to phase transformation problems. Numeri- 
cal results and analytical solutions were found to be in good agreement. 

A second aspect of this research dealt with the simulation of 
surface boundary temperatures by employing meteorological parameters 
that are readily available for most geographic locations of interest. 

It must be emphasized that the theoretical and experimental results of 
numerous investigators were relied upon to evaluate the individual 
hydrologic energy exchange terms discussed in Chapter 3. Simulated 
energy terms and surface boundary temperatures were the results of 
simplified approximations to complex thermal processes. The numerical 
simulations of surface energy exchange under Fairbanks summer and winter 
climatic conditions assisted in assessing the simplifying assumptions 
employed. Of particular interest were the small discrepancies between 
simulated and measured soil column temperatures at the end of the 10 
day Fairbanks winter simulation, considering that ambient air tempera- 
ture fluctuated from -39°F to 16°F. 

The final phase of the work applied the thermal model to an 
engineering analysis of silt and gravel offshore drill sites. Based 
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on the thermal simulation results, the following conclusions are made 
subject to Barter Island, Alaska, climatological conditions: 

1. Under natural freezing conditions with negligible snow cover, 
equilibrium frost penetration depth for saturated silt with 

a 4 foot gravel overburden is approximately 27 feet and for 
saturated gravel, 48 feet. 

2. Summer thaw of saturated frozen gravel at a porosity of .45 
can be expected to be 3 feet and at a porosity of .30 can 
be expected to be 5 feet. 

3. Dredged island drill sites are thermally feasible in shallow 
water depths. 

Expanding on the above conclusions, the depth of water a dredged 
island could be built in would depend on a number of variables such as 
the magnitude of lateral sea ice forces and the manner in which the 
island is dredged. A possible method of constructing a dredged island 
is to place the silt or gravel in an expendable metal torus perforated 
to relieve excess pore water pressures as the island freezes. Use of 
a torus would prevent loss of material due to wave erosion and hasten 
freezing in contrast to dredged soils deposited at a small natural 
angle of repose. Since both silt and gravel freeze faster than seawater, 
the metal torus would not be destroyed unles the compressive strength 
of the frozen materials was exceeded or the frictional resistance of 
the entire structure was less than the maximum expected lateral force 
due to sea ice movement. 

Lateral resistance to sea ice could be enhanced by driving piles 
at the center of the torus and by using sea bed foundation design methods 
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currently employed on gravity structures in the North Sea. The 
1 addition of an inclined rubble breakwater around the island would 
permit rafting of sea ice and reduce the horizontal component of ice 
forces. Furthermore, where drill sites are built in extremely shallow 
water, it becomes possible for seafloor sediments beneath the sites to 
freeze and increase resistance to lateral movement. 

5.2 Recommendations for Future Research 

Thermal analysis of permafrost island drill sites is no longer a 
problem of academic interest alone. Several drill sites have been 
constructed and occupied by Imperial Oil Limited of Canada in very 
shallow portions of the Mackenzie Delta (Riley, 1974), and plans are 
being developed to construct a dredged island drill site in 25 feet of 
water during the summer of 1974. As arctic offshore exploration 
progresses into deeper water, much larger ice forces will be encountered 
relative to the nearshore areas where ice mobility is somewhat limited. 

It will soon become necessary to develop methods of constructing in a 
relatively short time span massive structures capable of withstanding 
substantial ice sheet pressures. An economically feasible structure 
worthy of design consideration would be a grounded ice or ice-aggregate 
island. 

Before any large scale prototype ice island is built, a considerable 
amount of engineering research effort will be directed toward solving 
the thermal problem of ice accretion. An optimum water thickness will 
be sought which will be easy to apply, will minimize thermal degradation 
of ice already produced, and when repeatedly applied, will yield the 
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desired thickness of ice within the available time constraints. Of 
secondary interest will be the problem of determining the volume of 
seafloor sediments capable of freezing beneath the base of an ice 
island. This is of value in computing resistance to horizontal forces 
and estimating the size of the island. 

It is recommended that the finite element method be considered as 
a research tool for these problems. With minor modifications to the 
numerical model formulated in this report, useful design solutions to 
the ice accretion problem could be generated. No modifications are 
necessary to solve the seafloor sediment freezing problem. In addition 
to the thermal problems, it is also recommended that laboratory tests 
be conducted to evaluate freezing front boundary shear strengths of 
Beaufort Sea sediments at probable future offshore drill sites. 
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Variational Solution to the Heat Conduction Equation 

It is stated in Chapter 2 that, upon minimization, the generalized 
functional 



x = 7 / ’£ {k x c !t ) 2 + S^ 2 + 2Cpc H ■ G ^ }dxd x + 
♦ ±J>0f> 2 - 2<t> 0O <f>)ds + a/ Ci£ B 4> 5 - e A <J>i<t>)ds 



(A.l) 



where c is that portion of the boundary surface where values of <p are 
not imposed. 

satisfies the heat conduction equation and the associated boundary 
conditions of radiative, convective, and specified heat flux. The 
transition from equation (A.l) to the heat conduction equation will be 
demonstrated. 

Assume at a given instant in time 3^ is invariant and let the 

c3t 

first term in equation (A.l) be represented by any arbitrary functional 
of the form 

/£f (x,y,cf>,4> x ,4> y )dxdy (A. 2) 

Taking an arbitrary variation of f and its derivatives, produces 



6x ~ / £ [c lf 6<,>) + c H^ s< * , x ) + + (E q6< * ,ds 

x v (A. 3) 

+ jCh(<j)5<j) - 4> 0O 6c())ds + o£ (e B <t> 4 <5 <P - ewfwStfOds 

Because the 6 and operators commute and a minimum of the first 

dX 

variation is desired, the previous equation becomes 

■ wW* * & S§«« * % y 4 mmxiy 

(A.4) 

+ + h(<t><54> - QJty) + a(e B <p 4 6<p - ^is^Jds = 0 
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The last two terms in the first integral of equation CA.4) may be 
integrated by applying Green's theorem of integral calculus such that 



^ + - ''I'*'*®? * w^ ]ixiy 



vy 

* * ly % ,ds 



CA.5) 



where lx and ly are the direction cosines in the x and y directions 
respectively 

Upon substitution equation CA.5) into equation (A-4) we have 



& - &%? ]axiy 



CA.6) 



+ ^64) [ix||— + ly||- + q + h(4> - <J>co) + cr Ce B 4> 4 - £co<td)]ds = 0 
x r y 

Since <5<}> may be any arbitrary variation in R and <p is not imposed on c, 
application of the fundamental lemma of variational calculus requires 
that in R 



9f _ 3.9f _ 3 f 3f 

” 9x^3<J> x ” 9y*-9<J>y^ 



= 0 



(A. 7) 



and along c 

lx||- + lyff” + q + h C4> - 4>oo) + 0 (e B <f) 4 - e a,<t>i) = 0 (A. 8) 

x y 

By setting f equal to the first term of equation (A.l), it can be 
seen that we have an equivalent solution to the heat conduction equation 
if <5x is indeed a minimum. 
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User Instructions 

This finite element program solves one and two-dimensional 
transient heat conduction problems involving isothermal phase trans- 
formation. Boundary temperatures may be specified or internally genera- 
ted at the air-surface interface by providing specified meteorological 
parameters. One hundred and fifty elements and 100 nodes are permitted 
with the restriction that the numerical difference between any two 
numbered nodes common to an element not exceed 19; i.e. , a maximum 
system bandwidth of 20 is permitted. All nodes must be consecutively 
numbered beginning with one and all elements must remain within the 
first quadrant of the Cartesian coordinate system. If one-dimensional 
elements are used, they are required to be vertically arranged. 

Energy simulation across a two-dimensional surface is accomplished 
by the use of one-dimensional elements as discussed in Chapter 3. The 
one-dimensional surface nodes and elements must be consecutively 
numbered from left to right beginning with one. The remaining two- 
dimensional elements and their respective nodes may be numbered in any 
desirable sequence. Do not specify imposed boundary temperatures on any 
of the one-dimensional nodes used for surface energy simulation. 
Equivalent temperatures are computed and imposed within the program 
itself. During data input assign a zero volumetric heat capacity and 
any fictitious conductivity to the one-dimensional elements. New 
conductivities are computed after meteorological data has been read in. 

Thermal units of BTU/FT/DAY and °F must be employed if the surface 
energy simulation mode of the program is selected. Any system of 
consistent units may be used for imposed boundary value problems with 
the exception that time be expressed in days. 
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Data Input 

A. FIRST CONTROL CARD (613, 5F10.2) 



CC. 1-3 


1 for one-dimensional finite element simulation; 

2 for two-dimensional finite element simulation; 

3 for coupled one and two-dimensional elements for 
surface energy simulation. 


CC. 4-6 


Number of elements. 


CC. 7-9 


Number of specified temperature boundary nodes. 


CC. 10-12 


Number of nodes not held at a specified temperature 
less number of nodes used for surface energy simulation 


CC. 13-15 


Total number of nodes. 


CC. 16-18 


Number of one-dimensional elements used for surface 
energy simulation. 


CC. 19-28 


Latent heat of transformation per unit weight of 
substance that is subject to phase transformation. 


CC. 29-38 


Temperature at which isothermal phase transformation 
occurs. 


CC. 39-48 


Time step size in DAYS or decimal fractions thereof. 


CC. 49-58 


Desired elapsed time between nodal temperature print- 
outs, expressed in days. 


CC. 59-68 


Total simulation duration in days. 



B. SECOND CONTROL CARD (IS) 



CC. 1-S 



Number of elapsed time segments before new boundary 
or meteorological data is to be read in. 
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C. THIRD CONTROL CARD (1F4.1); (Delete if surface energy simulation 
t ' not desired) 



CC. 1-4 


Average snow depth present. If no snow cover. is being 
simulated, enter a zero. 



D. FOURTH CONTROL CARD (1F5.2); (Delete if surface energy simulation 
not desired) 



CC. 1-5 


Enter 1 if it is desired to reform system matrices in 
every time increment due to large fluctuations in 
meteorological parameters; otherwise enter a zero. 



E. FIRST DATA CARD (8F7.2) 



CC. 1-7 


List initial temperature of nodes, excluding specified 
temperature boundary nodes, in node number sequence. 


CC. 8-56 


Repeat above procedure. Use additional cards as 




necessary. 



F. SECOND DATA CARD (2F10.2) 



CC. 1-10 


X coordinate of node. 


CC. 11-20 


Y coordinate of node. 



CARDS MUST BE IN NODE NUMBER SEQUENCE! 
G. THIRD DATA CARD (313, 8F5.1) 



CC. 1-3 


First node number of element. 


CC. 4-6 


Second node number of element. 


CC. 7-9 


Third node number of element. Enter a zero if the 
element is one-dimensional. 

Node numbers of two-dimensional elements must be read 



in a counter clockwise direction around the element. 
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CC. 10-14 


Existing thermal conductivity of element. 


CC. 15-19 


Thawed thermal conductivity of element. 


CC. 20-24 


Frozen thermal conductivity of element. 


CC. 25-29 


Existing volumetric heat capacity of element. 


CC. 30-34 


Frozen volumetric heat capacity of element. 


CC. 35-39 


Thawed volumetric heat capacity of element. 


CC. 40-44 


Per cent moisture content of element. (Per cent of 
element unit weight). 


CC. 45-49 


Element unit weight. 



ALL CARDS MUST BE IN ELEMENT NUMBER SEQUENCE! 
H. FOURTH DATA CARD (7(I3,F7.2)) 



CC. 1-3 


Node number. 


CC. 4-10 


Specified boundary temperature. 


CC. 11-70 


Repeat above sequence. 



DATA MUST BE PROVIDED IN NODE SEQUENCE! 

I. FIFTH DATA CARD (I6,2X,7F8.2); (Delete if surface energy simulation 
not desired) 



CC. 1-6 


Date (day, month, year), i.e., 071271. 


CC. 9-16 


Ambient air temperature (°F) . 


CC. 17-24 


Wind velocity in miles per hour. 


CC. 25-32 


Cloud cover in tenths. 


CC. 33-40 


Cloud base altitude in 1000 feet. Enter 34 if 
ceiling is unlimited. 


CC. 41-48 


Relative humidity in per cent. 
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CC. 49-56 Clear day solar heat flux (BTU/ft^-day) . 

CC. 57-64 Outgoing evaporative heat flux (BTU/ft^-day) . 

J. REPEAT H and I based on the number of times data is to be updated 
or boundary conditions changed. 

Computer Program Listing 

Complete computer program listing and example output are provided 
at the end of this appendix. Printed energy balance terms of evaporative 
heat flux, net long wave radiation exchange, and convective heat flux 
are considered as heat losses from the air-surface interface when 



positive. 
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INTEGER GATE 

7000 FORMAT (6I3.5F10.2) 

7025 FORMAT (15) 

7050 FORMAT ( IOX « • OATA ERROR- UNSPECIFIED NODES* BOUNDARY NODES DO 
1 NOT EQUAL TOTAL NUMBER OF NOOES' I - 
7075 FORMAT (1FA.1) 

FORMAT ( 8F7.2) 

FORMAT (2F10.2) ' 

F0RMAT(3I3»8F5.1) 

FORMAT ( 7 ( I 3 * F7 .2 ) ) 

FORMAT ( 16 *2X*7F8.2) 

FORMAT ( 10X» 38HDAT A ERROR-SYSTEM 8 AND 



7100 
7150 
7200 
7250 
'7275 
73 0 0 
7350 
7*00 
7*50 
7500 
' 7510 
7520 
7525 
7530 
7550 
7600 



WIDTH TOO LARGE) 
FORMAT( 10 X » 22HNEG ATI VE AREA ELEMENT, 15) 

FORMAT ( IOX , AHNOQE »2X, I5»2X»13HFR0Z£N ON DAY,2X,F7.2>' 
FORMAT ( 10X,*HN0DE,2X,I5,2X,13HTHAWED ON DAY,2X,F7.2) 
FORMAT( IOX, ’COMPUTED NODAL TEMPERATURES ON DAY* , 2X.F7.2// ) 



FORMAT (7(2X, ’NODE* 

FORMAT ( 2X » 7 ( • 

FORMAT ( / / ,2X , 7 ( 
FORMAT ( 2X , 7 I » 



3X i 



TEMP.’) ) 
»)) 



)) 



),//) 



SPECIFIED TIME SPAN HAS BEEN 



F0RMAT(7( 2X, I*,3X,F6.2) ) 

FORMAT (IOX, 'SIMULATION OVER 
"1 COMPLETED* ) 

8000 FORMAT ( 10X.3X, ’ SWVP’ ,3X , 3X, • QSOL ’ ,3X , 2X , ’ EVFLX* ,3X, 2X, ’QLRAD* , 
13X.3X, »QCONV» ,3X) 

'8050 FORMAT!// , IOX, ’CGMPUTED SATURATION VAPOR PRESSURE AND ENERGY 
1BALANCE TERMS CN DAY’,F7.2) 

8100 FORMAT (10X,5FI0.2) 

8150 FORMAT ( IOX ,’ SNOW COVER HAS "MELTED SNOW COVER SIMULATION TERMINATE 
ID ON DAY* ,F7. 2,’ .',/, IOX, ’COMPUTED NODAL TEMPERATURES ARE ’ ) 

8175 FORMAT ( IOX, ’SIMULATION CALENCAR OATE IS(YRMOOY) ’,2X,I6,//) 

8200 FORMAT ( IOX, » SOLUT ION HAS BECCME UNSTABLE ENERGY SIMULATION HAS 
1BEEN TERMINATED’ ) 

8225 FORMAT ( IOX, • COMPUTED EQUIVALENT TEMPERATURE* , 1F6. 2, 13X ,’ AMB I ENT 
1AIR TEMPERATURE’ , 1F6. 2) 

82*0 FORMAT( // , IOX, ’THE FOLLOWING NOOES ARE PARTIALLY 

8250 FORMAT( IOX, ’NODE* ,5X, ’PERCENT OF ASSIGNED VOLUME 
1NED VOLUME’) 

827 5 FORMAT ( IOX, I* »19X»1F10.2»31X» 1F5.2) 



FROZEN •) 

FROZEN’ ,5X, ’ ASSIG 



8300 FORMAT( IOX, ’DATA ERROR 
1B0UNDARY CGNDITICN MUST 
8325 FORMAT!//, IOX, ’NO NUDES 

1 //) 

8350 FORMAT ( 1F5.2) 

COMMON /8LK 1/ R( 100) 

* COMMON /BLK2/ S( 100,20) 
COMMON /BLK3/ P( 100,20) " 
COMMON /8LK4/ SK 100,20) 
COMMON /BLK5/ NOD( 150,3) 
COMMON /BLK6/ X(100) 
COMMON /RLK7/ Y( 100) 
COMMUN /BLK8/ TK(150) 
COMMON /BLK9/ R0WC(150) 
COMMON / BLK10/ W( 100,20) 



NO BOUNDARY NODE SPECIFIED. AT LEAST ONE 
BE READ IN AS DATA’ ) 

UNDERGOING PHASE TRANSITION AT THIS TIME’ 



B-7 



COMMON /BLK 1 1 / Z(ICO) _ _ 

COMMON /BLK12/ NBC(50) 

COMMON /BLK13/ QC(50) 

DIMENSION RWCF ( LOO) , RWCT ( iOO) ,NKNC IOO) , TO ( 1 00 ) , G ( 100 > 
DIMENSION TXT ( 153) , TKF i I 50) , RCWC T ( 15 0 ) , ROWC F ( 1 50 ) , CMNGD ( 150 ) , 
IF LU I D ( 150 ) tFNDO ( 150) ,PCWTR( 150) ,UNWAT(150) 

DIMENSION VCL(IOO) ,FZHtT(100) ,NFREZ(100) »WATE( 100) 

DIMENSION A ( 3 , 3 ) , PEL ( 3 , 3 ) , XLCC ( 3 ) , YLOC( 3 ) , SEL ( 3,3 ) 

READ! 1,7000) KODE , NEMS , L BC , NK , MM , NSURF , HE AT , TF AZ , DT ,0 AY , END 
READ(1,7025) NTSEG 

COUNT=CAY/DT . ~ ' 

C0UNTI=0. 

TIME = 0. 

TCAY= 0. ' *“■■■ * ~~ 

ITSEG = 0. 

SNOW - 0. ' _ _ 

” OPT = 0. 

IF ( NSURF ) 200,200,100 
100 READ ( 1,7075) SNOW 
READ! 1,8350) OPT 
200 CONTINUE 

I F (MM- ( NK+L8C ) .EQ. 0) GO TO 500 
WR I TE { 3 , 7C50 ) 

GO TO 3900 
500 CONTINUE 



C READ INITAL NODAL TEMPERATURES 

READ (T77 ICO) ~C~ TO ( I ) “ 7 T= 1,NK ) 



00 1000 1=1, NK 
VOL ( I) - 0* 

WATE(I) = 0. 

RWCF { I) = 0. 

RWCT ( I ) = 0. 

IFITOII) *.GT. TFA'ZT’GO T0“~ 800 
NFREZ(I) =1 
GO TO 1C00 
800 MFREZm = 2 
1000 CONTINUE 



C— READ FINITE ELEMENT MESH COORDINATES AND ELEMENT THERMAL 
C PARAMETERS 



R E AO (1,7150) C X C I),Y( I) 
READ ( 1,7200) ( (NOD(I,J 

1R0WCFI I) , RGWC T ( l ) ,PCWTR 
DD 1005 I = 1, NEMS 
FNQD ( I ) = 0 • 

1005 CMNOD(I) = 0* 

. 1006 IF ( LBC ) 4000 , 4000 , 1007 
1007 NSERT = NSURF ♦ 1 

_ IF ( NSERT .GT. LBC) GO 
READ( 1,7250, END =3800) 
1113 IF C NSURF) 10C8 , ICO 3 , 1 



, 1=1, MM) 

) , J= 1 » 3 ) ,TK( I 
( I) , U N W A T ( I) , 



TO 1113 
( NBC ( J ) , BC ( J ) 
111 



) , TKT ( I ) , TKF ( I ) , ROWC { I), 
1=1, NEMS) 



, J = NSERT , LBC ) 



B-8 



f 



.1111 DO 1112 1=1, NSURF 

1112 NBC ( l ) = I 
C 

2 ______ C DETERMIN E U NKNOWN NODE _P01_NTS 

1008 IF (TIME) 3900* *1009,1355 

1009 I = 1 

' NKNC I ) =1 

1010 DO 1050 J = 1 , L8C 
LL=NBC ( J ) 

IF ( LL-NKN(I) VnE.’OI GO TO' 1050 
NKN ( I ) = NKN ( I ) + 1 
1050 CONTINUE 
' 1060 1=1+1 

IF ( I .GT. NK) GO TO 1070 
NKN ( I ) =NKN( I— 1 ) + 1 
" GO TO 1010' 

1070 CONTINUE 



C— —COMPUTE THE NUMBER OF UNSPECIFIED NODES COMMON TO EACH ELEMENT 



IF ( NSURF .LE. 0 ) GO TO 1300 
DO 1225 1 = 1, NSURF ~ ” 

1225 CMNOD ( I) = 1. 

NSERT=NSURF+ 1 

GO TO 1325 

1300 NSER T= 1 

1325 DO 135C J=NSERT , NEMS 
... DQ l35Q K = 1>3 

L= 0 

001350 M=1,LBC 

IF (MBC ( M ) — NOD ( J , K ) .EQ. 0) GO TO 1350 
IF ( NOD ( J , K ) .EQ. 0 ) GO TO 1350 
L = L + 1 

IF < C~.LTrLBC) GO TO 1350 
CMNOD ( J ) = CMNOD ( J ) + 1. 

1350 CONTINUE 
1355 CONTINUE 

IF (NSURF) 1300,1800,1360 

1360 READ! 1 ,7275, END=3800) DATE , TAMe, VE L, CLOUD, ALT, RH , SRFLX ,6VFLX 
"IF (VEL .GT. 1.) GO TO 1380 
VEL = 1. 

C -ASSIGN APPROPRIATE THERMAL AND ROUGHNESS PARAMETERS' 

1300 IF ( SNOW) 1385,1385,1395 
1385 EMISS = .90 

• ALBDO = .15 
GO TO 1430 
1395 EMISS = .98 

IF ( TAMB - 32.0) 1400,1405,1405 

1400 ALBDO = .85 _ 

' GO TO 1420 " 

1405 IF ( TO AY ) 1410,1410,1415 



o o o 



1410 TOAY = TIME 

1415 AGE = ABSITIME - TOAY) 

AIBOO = .62*EXP(-.0215*AGE) ♦ .23*EXP(-.49*AGE) 
1420 CONTINUE 



COMPUTE CONVECTIVE HEAT TRANSFER COEFFICIENT. ! 

L m , ^ _ J 

1430 IF I ABSITAVE'- TAMB) .LT. 1.) GO TO 1440 

IF (TAVE - TAMB) 1435,1440,1445 ; 

1435 IF ISNOW) 1436,1436,1437 

1436 CNVCT = 49.6*VEL*.38 ' "7 ’ 

GO TO 1800 

1437 CNVCT = 25.*VEL* .38 

GO TO 1800 " " 

1440 IF ( SNOW) 1441,1441,1442 

1441 CNVCT = 86.6*VEL*.30 

GO TO 1800 ' ‘ [ 

1442 CNVCT = 60.*VEL*.30 j 

GO TO 1800 

1445 IF I SNOW) 1446,1446,1447' - [ 

1446 CNVCT = 120.*VEL*.22 1 

GO TO 1800 ' | 

'1447 CNVCT = 97.*VEL*.22 7 

GO TO 1800 

1800 REDO * 0. 

IF (OPT .LT. 1. ) GO TO 1810 
REDO = 1. 

1810 CONTINUE 
NROW = MM 

IF (NSURF) 1803,1803,1801 ' v ;’ 

1801 SUM = C. 

00 1802 I - 1, NSURF ' " 

1802 SUM = SUM *■ TO ( I ) 1 

TAVE = SUM /FLOAT ( NSURF ) 

IF ( ABS ( TAVE -TAMB) .LT. 100.) GO TO' 1803 ' 

WRITE (3,8200) 

GO TO 3900 

1803 IF (TIME) 1805,1805,2020' 

1805 N0UM=1 

C— -COMPUTE SYSTEM - BAN0 _ W"I0TH 7 

NC0L=1 

IF ( KOOE .NE'. 1 ) GO TO 1900 

I F ( KODE .EO. 1) NC0L=2 
GO TO 2010 
1900 CONTINUE 

IFIKOOE .EQ. 2) GO TO 1950 
NDUM= NSURF +• 1 
1950 00 2000 N-NOUM,NEMS 
00 2000 1=1,3 

00 20C0 J= 1 » 3 

NN = NGO(N,I) - NOO(N,J) ♦ 1 
2000 IF ( NCOL - NN .LT. 0 ) MCOL=NN 
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2010 CONTINUE ; . 

IF ( NCOL .LE. 20 I GO TO 2020 
WRITE (3,73001 

GO TO 3 90 0 

C ZERO SYSTEM ARRAYS 

C 

2020 DO 2040 1*1, MM 
R ( I ) = 0. 

00 2040 J=1, NC0L ■ 

ST( I , J) = 0. 

S(I,J) = 0. 

P ( I , J I = 0. 

2040" W ( I , J ) = 0. 

IF {NSURF I 2045,2045,2041 

C 

C — - — COMPUTE INCIDENT SOLAR FLUX AND' LONG WAVE ENERGY EXCHANGE 

2041 QSOLM 1.— ALB DO I* ( l.-( . 82-. 024*AL T) *CL0UD/10. ) *SRFLX 
TABS = (TAMB-32. )*<5./9. ) ♦ 273.16 

TS=( TAVE— 32. ) * ( 5 . /9. ) ♦ 273.16 
VALUE = 373.16/TABS 

SVPL = — 7 .9029 8* ( VALUE —1.1 ♦ 5 . 02808*ALOGl0( VALUE I 
l-( 1.3816E-7) *< 10.**{ IL.344M 1 .-( 1 . /VALUE ) ) )-l. I 
2+ ( 8 . 13 2 BE— 3 )*( 10.**( -3.45149* (VALUE -l.ll -l.J 
3 ♦ 3.0C5715 ' 

SWVP= 10.**(SVPL I 
VP=SWVP*RH/100. 

— "TRANS= 1.0 — . 024* ALT 

ATMEM = .74 + .0125*CL0UD ♦ .0049 *VP 

QLRAD= 1 1 . — TRANS*CL0UD/10 .) *(4.389E-7)*(EMISS*(T S**4. J-ATMEM* 
1 ( TABS**4. ) ) 

_C- C ONVERT RA DI ATION F LUXES TO AN EQUIVALENT TEMPERATURE 

TEQ = (QSOL - EVFLX - QLRAD) /( . 336*CNVCT) 

EQTEM = TAMB ♦ TEQ_ 

2042 DO 2043 I=1,NSURF 
BC (II = EQTEM 

2043 TK ( I ) = .336*CNVC T 

C CONSTRUCT THE SYSTEM MATRICES 

2045 D0~2600~ M="l"i NE MS 

CALL CENTER ( M, X6AR, Y3AR .NSURF, KODE, DELTY, AM) 

IF ( KCDE .EQ. 1) GO TO 2550 
IF ( M~ .LE . NSURF) GU TO 2550 
CALL XIETA(M,XBAR, YBAR , XLOC , YLOC ) 

CALL AREA ( XLOC , YLOC , AM) 

IF ( AM .GT. 0.) GO TO 2550 

WRITE(3,7350) M 
GO TO 3900 
2550 CONTINUE 

2050 IF ( M .LE. NSURF I GO TO 2575 



) 
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IF ( TIME .GT. 0.) GO TO 2575 

00 2200 K = 1,3 
J * 0 

00 22CC L * 1,L8C 
NCOM '= NODIMVK) 

IF ( NBC ( L ) - NCOM .NE. 0) GO TO 2075 
GO TO 2200 

2075 IF I NCCM .EQ. 0 ) GO TO 2200 
J = J ♦ l 

IF (J .EQ. LBC ) GO TO 2100 
GO' "T0‘ 2200 

2100 IF ( KCDE .EQ. 1 ) GO TO 2110 
BETA = AM 

GO TO 2120 ' 

2110 BETA = ABS (DELTY)*AM 
2120 DO 2125 I = l.NK 

IF ( NKMI) — NCGM .EQ. 0 T GO TO 2130 " 

2125 CONTINUE 

2130 FLUID ( M ) = BETA*PCWTR(M)*UNWAT(M)*(.01) 

VOUI) = VOL ( I ) + 3ETA/CMN0DI M) ' 

VlATE(I) = WATE(I) ♦ FLUID! M) /CMNOD(M) 

RWCT(I) = RWCT(I) * ( ROWCT { M ) *BE T A ) /CMNODI M ) 

RWCF(I) = RWCF(I) «- ( ROWCF ( M ) *BET A ) /CMNOD (Ml 
IF( TO(I) .GE. TFAZ) GO TO 2200 
FNOD(M) = FNOD(M) ♦ l. 

'2200 CONTINUE' 

2575 CONTINUE 

CALL S W ALS(M»SEL,AM,A,NSURF, KCDE , DELT Y) 

CALL SMALP ( M, AM ♦ PEL, NSURF , KODE, D ELTY ) 

CALL ADDIN ( SEL, PEL,M, KODE ) 

2600 CONTINUE 

IF'(TIME)~3900» 2610,2650 
2610 CONTINUE 

DO 2625 1=1, NK 

FZHET ( 11"= WATECI ) *HEA‘T : 

RWCF(I) = RWCF ( I ) /VOL ( I ) 

2625 RWCT(I) = RWCTm/VOLII) 

2650 CONTINUE "" ' 

C 

C INSERT SPECIFIED BOUNDARY TEMPERATURES AND MODIFY THE S,P, AND 

C— ---R MATRICES ACCORDINGLY 

C_ALL BNDINJN^ROW^NCOL.LBC) __ 

C BEGIN TIME DOMAIN SOLUTION (CRANK NICOLSON) 

C 

00 '2700" 1 = 1VN.R0 W 

DO 2700 J=1,NC0L 

W(I,J)= S(I,J) (2./DT)*P(I,J) 

2700 P(I,J) = (2./DT )*( P( I , J ) ) - S ( I , J ) 

C^^“PRETRIANGULARIZE_THE W MATfUX 



CALL PRESCLINROW, NCOL ) 
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2800 CONTINUE 

CALL CCM8(T0»NR0W»NC0L) 

DO 2900 I® 1 » NRGW 
2900 G ( I) = 2. *R( I) + III) 

"CALL FINSOL I G »NROW» NCOL ) 

TIME * TIME ♦ DT 
COUNT I = COUNT I * l 

ITS6G IT SEG ' ♦ I " * 

1*1 

GO TO 308 A 

3081 G(I) = TFAZ 

3082 TO (I) = G(I> 

1 = 1 ♦ I 

30 8 A CONTINUE ' 

IF ( I .GT. NROW ) GO TO 3689 

DIFF = G ( I ) - TFAZ 

IF ( DIFF) 3085 »3C 81 i 3487 

C 

C- — — SIMULAT E FR EEZING PROCESS 

3085 IF ( NFREZ(I) .EQ. I ) GO TO 3082 
DELQ = RWCTI I)*VUL( I)*DIFF 
FZHETII) ® FZHETII) ♦ DELQ 
IF { FZHETII) .GT. 0.) GO TO 3081 
NFREZ ( I) = 1 

G ( I ) = TFAZ FZHETII) / I RWCFI I )*VOL( I ) ) 

WR I TE ( 3 t 7400 ) NKN ( I ) * TIME 
FZHET (I) = WATE I I )*HEAT 

GO TO 3AOO * 

C 

C SIMULATE THAWING PROCESS 

3A87 IF ( NFREZII) .NE. I) GO TO 3C82 
3500 DELTQ = RWCFI I)*VOL( I)*DlFF*(-I.O) 
i FZHETI I) = FZHETII) + DELTQ 
IF (FZHETII) .GT. 0. ) GO TO 3081 
NFREZII) = 2 

G( I ) = TFAZ - FZHETII) / I RWCTI I )*VOL I I ) ) " 

WRl TE I 3 » 7 A50 ) NKN( I )» TIME 
FZHETII) = WATEI I)«HEAT 
3A00 NOON' = NKN I I ) ’ 

DO 3AI5 J= 1 » NEMS 
DO 3AI5 K=1 i 3 
"NCHOS® NODIJ.K) 

IF I NCHOS .EQ. NODN) GO TO 3AI0 
GO TO 3AI5 

34 1 0 - 1 F I ~ j.LE~NSURF)~GO‘'TO“ 3A15 

IF I NFREZII) .EQ. 2 ) GO TO 3412 
FNOD(J) = FNODIJ) 1. 

GO TO 3413 

3412 FNODIJ) = FNODIJ) - I. 

C-— -MODI FY - ELEMENT' THERMAL PARAMETERS WHEN NODAL PHASE CHANGE OCCURS - 
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3413 

3415 



TK(J) = ( 
ROWC(J) = 
1CHNOD ( J ) 
CONTINUE 
' REOO = 1.0 
GO TO 3082 
IF (NSURF) 



TKF ( J ) *F NOD ( J ) * ( 
( ROWCF(J)*FNQD(J) 



CHNCO ( J )-FNOD ( J ) ) *TKT ( J ) ) /CMNOD( J ) 

+ l CMNOD( J ) -FNODi J ) )*ROWCT(J) )/ 



3694, 36S4, 3690 



3689 

C COMPUTE AVERAGE SURFACE TEMPERATURE 

~~ 3690 



3691 



3692 

3693 



3694 



3695 



"3700 

3710 

3720 



3800 



3801 



SUM = 0. 

00 3691 I =1 , NSURF 
SUM = SUM + TOC I I 
TAVE = SUM/FLOAT (NSURF) 

QCONV = CNVCT*.336*(TAVE - TAM8) 

I F ( SNOW ) 3694,3694,3692 

IF ( TAVE - 32.) 3694,3693,3693 



TIME 



( NKN(I), TOC I ) , I = 1 ,NK) 



( NKNC I ) , TOC I ) , 1=1, NK) 



SNOW = 0. 

WRITE(3,8150) 

WRI TE ( 3, 7525 ) 

KR I TE( 3 , 75 50 ) 

WR I TE ( 3 , 7520 ) 

GO TO 3900 ‘ 

IF (TIME .GE. END) GO TO 38G0 
IF (COUNTI .LT. COUNT) GO TO 3695 
COUNTI = 0/ 

WRITE ( 3,7500) TIME 
WR I TE ( 3,7520) 

WRI TEC 3, 1 5 10) 

WRITE(3,7520) 

WRITEC3 ,7550) 

WR I TE ( 3 , 7 520 ) 

IF (NSURF .LT. 

WR I TE ( 3 , 8050 ) 

WRITE(3,8175) 

WR I TE( 3 , 8225 ) 

WR I TE ( 3,7520) 

WR ITE ( 3 ,8000 ) 

WRITE(3,7520) 

WR I TE ( 3 , 8100 ) 

WRITE (3,7530) 

IF ( I TSEG .LT. 

ITSEG = 0 

IF (ITSEG)' 3710,3710,3720 
GO TO 1006 

IF (RECO .GE. 1.) GO TO 1380 
GO TO 2800 ‘ 

WRITE (3,7600) 



► l ) GO 
TIME 
OATE 



TO 3695 



EQTEM, TAM8 



SWVP ,QSOL , EVFLX ,QLR AO, QCONV 



NTSEG) GO TO 3700 



WR I TE ( 3 , 7 500 ) 
IF (NSURF. LT. 
WRITE(3,8175) 
WRITE(3,7520) 
WR I TE ( 3 , 75 10 ) 
WRITE(3,7520) 



TIME 
1) GO 
DATE 



TO 3801 



t 
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— WR I TE ( 3 » 7 5 50 ) C NKN( I } , TQ( I ) , 1=1, NK) 

WRITE (3,7520) 

MNTR= 0 
CO 3810 1 = 1, NK 

IF ( TO ( I ) - TFAZ ) '3810, 3804, '3810 ' 
3804 MNTR = MNTR ♦ 1 

IF I MNTR »NE _GO TO 38C7 
WR I TE I 3 ,8 240 ) 

WRITEI3.7525) 

WRITE (3, *8250) 

WR I TE ( 3 , 7 520 ) 

3807 CONTINUE 

CK= FZHET ( I )/(WATE(I )*HE AT ) 

IF ( NFREZ(I) .EQ. 1) GO TO 3808 
PRCNT = ( 1. - CK 1*100. 

GO TO 3809 

'3808' PRCNT "= CK*'l00. 

3809 WR I TE ( 3 , 8275 ) NKN( I ) , PRCNT, VOLIII . 

3810 CONTINUE 

'IF (MNTR .GT~. 6) GO TO 3811 ' 

WRITEI3.8325) 

3811 WRITE(3,7530) 

3900 CONTINUE 

GO TO 5000 
4000 WR I TE ( 3 ,8300 ) 

5000 CONTINUE ' 

STOP 

END 



4 



Moor* tuiln*i« Form*, Inc. »* 



B-15 



* 



C 

C 



...SUBROUTINE CENTER ( M , X BAR ,YBAR , NSURF , KOD E ,D ELTY , AM ) _ 

THIS SUBROUTINE COMPUTES THE CENTER COORDINATE OF TRIANGULAR 
ELEMENT M. 

COMMON /BLK 5/ N0D(150»3) 

COMMON /BLK6/ X ( IOO ) 

COMMON /BLK 7/ Y(IOO) 

_I_= NOD I M, I) 

J = N0D(M,2T 
IF I KODE.EO. I) GO TO 2 
_ IF (KODE .EO. 2 ) GO TO 
IF ( M .LE. NSURF ) GO TO 

1 K = NOD ( M ,3 ) 

XBAR _= ( X(I) + X(J) + X ( K ) 

YBAR = '( YU) + YIJ) + Y ( K ) 

GO TO 9 

2 DELTY = YIJ) - Y ( I ) 

AM = 1 . 

IF ( KODE .NE. 3 ) 

AM =_ I X ( NSURF ) - 

'9 CONTINUE 
RETURN 
END 



1 

2 



) / 3. 
) /3. 



GO TO 9 
XU) ) / 



FLOAT! NSURF) 



( 



1 



( 
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SUBROUTINE XIETA { M ,XBAR , YB AR ,X LOC , YLOC > 

C THIS SUBROUTINE COMPUTES THE LOC A L_ COORDINATES OF 

C NODES ON TRIANGULAR ELEMENTS M. 

COMMON /RLK5/ N0D(150t3) 

COMMON /BLK 6/ X(IOO) ” 

COMMON /RLK7/ Y(IOO) 

DIMENSION XLOC ( 3 ) » YLOC( 3) 

00 100 L=1 » 3 ; 

1= NOD ( M» L ) 

XLOC ( L ) = X (I) -X BAR 
100 YLOC < L ) = Y ( I ) — YBAR 
RETURN ' ' 

END 



( 

L 






V. 
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SUBROUTINE area' (XLOC",YLOC,AMf 

THIS SUBROUTINE CALCULATES THE AREA OF TRIANGULAR 
ELEMENT M. 

DIMENSION XL0C(3)»YL0C(3)‘ 

AM= ( XLOC ( 2 ) *YLOC ( 3 ) — XLOC ( 3 ) *Y LOC ( 2 ) -XLOC ( 1) *YLOC(3) 

1 +XL0C(3)*YL0C( 1) +XLOC(l)*YLOC(2)-XLOC(2>*YLOC< 1) J/2. 
RETURN 
END 



(' 



( 



i. 






f 



( 



V 



( 
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SUBROUTINE SM AL S ( M , SEL , AM, A , NSURF , KOOE , OF LT Y) 

C THIS SUBROUTINE COMPUTES THE SMALL S MATRIX 

C FOR TRIANGULAR ELEMENT M. 

_ COMMON /BLK 5/ N0D(150,3) 

f COMMON /BLK6/ X(IOO) 

I COMMON /BLK7/ Y(IOO) 

COMMON /BLK8/ TK(150) 

. DIMENSION SEL (3 ,3) ,A( 3,3) 

I = NOD (M ,1 ) 

J =_ NOD (M* 2 ) • 

IF ( KOOE .EO. 1 ) GO TO 3 
IF ( KODE .EO. 2 ) GO TO 1 
I F _ ( M .LE. NSURF) GO TO 3 

1 K = NOD ( M * 3 ) 

A (1,1 >=X(J)*Y(K)-X(K)*Y( J) 

_A( 1,2>=X(K )*Y(I ) —X ( I )*Y(K) 

A ( 1 , 3 ) =X ( I ) *Y ( J ) -X ( J ) *Y ( I ) 

A ( 2 , 1 ) = Y ( J ) — Y ( K ) 

A ( 2 ,2 ) = Y (K ) — Y ( I ) 

A ( 2 » 3 ) = Y ( I ) — Y ( J ) 

A( 3,1) = X ( K ) - X( J) 

A< 3,2 ) = X ( I ) - X ( K > ^ 

A ( 3,3) = X(J ) - X ( I ) 

DO 2 J = 1 ,3 

DO 2 K = 1,3 

2 SEL ( J , K ) = A( 2 , J ) *A ("2 , K )'=*"( TK( M )/ ( AM*4 . ) ) 

1 + A ( 3 » J ) *A ( 3 * K ) * ( TK(M)/( AM*4. ) ) 

GO TO 5 

'3 ELMS = ( TK ( M )* AM ) /DELTY 

SEL ( 1 , 1 ) = ELMS*1. 

SEL ( 1 ,2 ) = ELMS* ( — 1. ) 

SEL(2 , 1 ) = ELMS*(-1.) 

SEL (2 ,2 ) = ELMS*1. 

DO A I = 1,3 

4 SEL ( I ,3 f = 0. 

5 CONTINUE 
RETURN 
END ' 
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SUBROUTINE SMALP (M, AM,PEL,NSURF ,KOOE,OFLTY_ ) 

C THIS SUBROUTINE CONSTRUCTS A SYMMETRICAL P ELEMENT 

n C SUBSYSTEM MATRIX FOR THE MTH TRIANGULAR ELEMENT 

COMMON /RLK9/ R0WCI150) 

DIMENSION PEL (3,3) 

IF ( KODE .EQ. 1 ) GO TO 2 

I F ( KODE .EO. 2 ) GO TO 1 

IF ( M .LE. NSURF) GOTO 2 

1 EL « ( ROWC C M ) AM ) / 12. 

PEL (1,1) = 2 .*EL 

PEL(1,2)=1.*EL 
PEL ( 1 » 3 ) =1 .*EL 

PEL (2,1 >=1.*EL 

PEL (2,2 ) =2 .*EL 
PEL (2 »3 ) =1 . *EL 
PEL (3 » 1 ) =1 .*EL 
PEL (3*2) =1 .*EL 
PEL (3»3)=2 .*EL 

GO TO 5 

2 ELP = ( AM*R0WC(M)*DELTY)/6. 

PEL ( 1 ,1) = ELP*2. 

PEL (1*2) = ELP*1. 

PEL (2*1) = ELP*1. 

PEL (2,2) = ELP*2. 

DO 3 1=1,3 

3 PEL ( I , 3 ) = 0. 

5 CONTINUE 

_ RETURN . . 

ENO 




(. 



SUBROUTINE ADD I N < SE L , PE L , M , KODEi 

COMMON /BIK2/ S ( 100» 20) 

COMMON /BLK5/ NOD( 150,3) 

COMMON /BLK 3/ P( 100,20) 

DIMENSION SEL ( 3 , 3 ) 

DIMENSION PEL (3,3) 

IF_( KODE .EO. i) GO TO 2 

IF ( KODE .EO. '2 ) GO TO 1 
IF ( M .LE. NSURF ) GO TO 2 
IDO 100 J = 1,3 
NR = NOD (M , J ) 

DO 100 K = 1,3 

NC = NOD (M ,K ) - NOD ( M , J ) + 1 

IF ( NC .LE. O') GO TO 100 

S ( NR , NC ) = S ( NR ,NC ) ♦ SEL (J,K) 

P ( NR , NC ) = P ( NR , NC ) + PEL(J,K) 

100 ' CONTI NUE 

GO TO 5 

2 DO 200 J=1 ,2 
’ "NR = NOD ( M , J ) 

DO 200 K = 1,2 

NC = NOD ( M, K ) - NOD(M, J) + 1 
"IF ( NC .LE. 0) GO TO 200 
S< NR ,NC ) = S ( NR ,NC ) + SEL(J,K) 

P ( NR ,NC ) a P ( NR , NC ) + PEL (J,K) 
200 CONTINUE ~ 

5 CONTINUE 
RETURN 
"END 
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SUBROUTINE BNDIN ( NROW , NCOL , LBC ) 

COMMON /BLK1/ R ( 100) 

□ COMMON /BLK2/ S( 100,20) 

COMMON /BLK3/ P(100,20) 

j COMMON /BLK12/ NBC(50) 

( I COMMON /BLK1 3/ BC(50) 

C — 

C- SET BOUNDARY CONDITIONS IN MATRICES 
C- 

:_D0 3230 1 = 1, LBC 

NN=NBC ( I ) 

«. NC=NCOL 

3200_NR= NN-NC+1 

I F ( NR .LT. 1 ) GO TO 3210 
R ( NR ) = R ( NR ) - S(NR,NC)-BC(I) 

_ 32 10_ NC = NC- 1 

IF(NC .GT. 1) GO TO 3200 " 

NC= 1 

NR=NN 

3220 NR=NR+ 1 

IF (NR .GT. NROW) GO TO 3230 

NC = NC+1_ 

• IF (NC .GT. NCOlTGO TO 3230 
R ( NR ) = R(NR)-S(NN,NC)*BC(I ) 

GO TO 3220 
3230 CONTINUE 
C- 

C-R6F0RM MATRICES AND ELIMINATE EQUATIONS 

C-AT 'BOUNDARY CONDITION NODES" 

C- 

DO 3300 1=1, LBC 

NN=NBC fl ) -1+1 
NR0W=NR0W-1 

I F (NN .GE. NROW) G0_ TO 32^50 

DO 3240 NR=NN , NROW 
l R ( NR ) =R (NR +1) 

DO 3240 NC = 1 , NCOL 

S(NR,NC)=S(NR+1,NC)'~ 

3240 P(NR,NC)=P(NR+1,NC) 

3250 NR=NN 
N=2 ‘ 

3260 NR=NR- 1 

; IF(NR. LT; 1 ) GO TO 3300 

N=N+1 > 

I F ( N .GT. NCOL) GO TO 3280 
DO 3270 J=N,NCOL 
P ( NR , J— 1 ) =P ( NR , J ) 

3270 S(NR,J-1)=S(NR, J) 

P (NR, NCOL ) =0. 

S(NR ,NCOL ) =0." ' ' 

GO TO 3260 

3280_S(NR,NC0L)=0. 

P (NR, NCOL >=0. 

3300 CONTINUE 
RETURN 
END 
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SUBROUTINE _PRESOL (NROW,NCOL) 



C-THIS SUBROUTINE PRE -TRI ANGULAR I ZES A SYMMETRIC 
C-MATRIC IN BAND FORM FOR SOLUTION BY THE GAUSSIAN 
C-ELIMINATION METHOD. FINAL SOLUTION IS BY SUBROUTINE FINSOL. 

C- 

COMMON /BLK4/ ST( 100,20) 

COMMON /BLK10 / W( 100,20) 

N=0 

100 N=N+ 1 

STIN,NCOL) =_W(N,1) 

C- 

C-REDUCE PIVOT EQUATION 

C— _ 

IF(N-NROW) '150,300,150 
150 DO 200 K=2 ,NCOL 

ST( N,K-1 )=W (N,K)_ 

200 W I N » K ) = WIN, K) /WIN, 1 ) 

c- 

C-REOUCE REMAINING EQUATIONS 

c- ' ' ‘ " " ~ ~ ~ 

DO 260 L=2,NC0L 
I =N+L— 1 

IFINROW-I) 260,240, 2 AO 
2 AO J=0 

00 2 50 K=L ,NCOL 

J = J+1 

250 HI I ,J)-W(I ,J)-ST(N,L-1>*W(N,K) 

260 CONTINUE 

GOTO 100 

300 RETURN 
END 



l 



i 
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SUBROUTINE _COMB ( Y , NROW »NCOL ) 

C- 

n c 

C- THIS SUBROUTINE MULTIPLIES SYMMETRIC MATRIX P 

T C-TO MATRIX Y AND STORES THE RESULT IN Z. 

C- __ 

COMMON /BLK3/ P ('106720') 

COMMON /BLKI 1/ Z( 100) 

DIMENS ION Y(100) 

DO 200 I = 1 » NRDW 
Z ( I ) =Y ( I )*P(I,1) 

_00 200 K = 2 1 NCOL 
L=I-K+1 

IF(L.LT.l) GO TO 100 
Zm-Z(I) + Y(L)*P(L»K) 

'100 N=I+K-1 ' 

IF(N.GT.NROW) GO TO 200 
Z(I)=Z(I) +Y(N)*P(I,K) 

200 CONTINUE 
RETURN 
ENO 
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SUBROUTINE FINSOLI SS,NROW,NCOLJ 



C-THIS SUBROUTINE SOLVES A SET OF LINEAR SIMULTANEOUS 
C-EOUATIONS WHOSE COEFFICIENT MATRIX HAS BEEN PR E-TR I ANGULARI ZED 
C-BY THE GAUSSIAN ELIMINATION METHOD. THE SYSTEM MATRIX 
C-IS IN BAND FORM AND IS SYMMETRICAL. SOLUTION IS PLACED 
C— IN THE LOAD MATRIX. 

- c — — : — — ~ 

C- 

COMMON /RLK4/ ST( 100,20) 

COMMON /BLK 1 0/ W( 100,20) " 

DIMENSION SS(IOO) 

c- 

C-REDUCE LOAD VARIABLE'S " " 

C- 

N=0 

100“N=N+I 

SS ( N ) = SS (N ) /S T ( N , NCOL ) 

IF (N-NROW) 110,200,110 
YlO CONTINUE 

DO 130 K=2 ♦ NCOL 
L=N+ K— 1 

IF(NROW-L) 130,120,120 
120 SS(L )=SS(L )-ST(N,K-l)*SS(N) 

130 CONTINUE 
GO 'TO 100 
C- 

C-BACK SUBSTITUTION 

c- 

200 N=NROW 
300 N=N-1 

IF ( N ) 350,500,350 
350 DO 400 K=2,NC0L 
L=N+ K— 1 

IF(NROW-L) 400, 3 7 O', 3^7 O' 

370 SS ( N ) = SS ( N ) - W ( N » K )*SS ( L ) 

400 CONTINUE 
~ GO "TO 300 
500 RETURN 
END 



( 
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